library(tidyverse)
library(fixest)

# See the README for instructions on how to access this data
dat <- read_csv("exposure_mvmt_timing_regressions_long_all.csv")

dat <- dat_in %>%
  mutate(
    fw_median_hh_inc_grp = cut_number(frnd_weighted_median_hh_inc, 100),
    fw_pop_density_grp = cut_number(frnd_weighted_pop_density, 100),
    fw_frac_urban_grp = cut_number(frnd_weighted_frac_urban, 100)) %>%
  mutate(
    fw_median_hh_inc_grp_X_month = paste0(as.numeric(fw_median_hh_inc_grp), "__", month),
    fw_pop_density_grp_X_month = paste0(as.numeric(fw_pop_density_grp), "__", month),
    fw_frac_urban_grp_X_month = paste0(as.numeric(fw_frac_urban_grp), "__", month)
  ) %>%
  mutate(
    zcta_interaction_group_X_month = paste0(home_zcta, "__", is_female, "__", has_college, "__",
                                            age_bracket, "__", has_iphone, "__", has_tablet, "__", month))


# This is a sanity regression that should simply match our baseline table run using reghdfe
ols_sanity <- fixest::feglm(d_share_no_mvmt ~ d_log_fwc_100k | as.factor(userid) +
                              fw_median_hh_inc_grp_X_month + fw_pop_density_grp_X_month + fw_frac_urban_grp_X_month +
                              zcta_interaction_group_X_month,
                            data=dat, fixef.rm = "both", family = gaussian())

# These are the new robustness regressions (Table A11)
ols_tiles_diff_no_userid <- fixest::feglm(pct_cg_avg_tiles ~ d_log_fwc_100k | fw_median_hh_inc_grp_X_month + fw_pop_density_grp_X_month + fw_frac_urban_grp_X_month +
                             zcta_interaction_group_X_month,
                           data=dat, fixef.rm = "both", family = gaussian())

ols_tiles_diff <- fixest::feglm(pct_cg_avg_tiles ~ d_log_fwc_100k | as.factor(userid) +
                             fw_median_hh_inc_grp_X_month + fw_pop_density_grp_X_month + fw_frac_urban_grp_X_month +
                             zcta_interaction_group_X_month,
                           data=dat, fixef.rm = "both", family = gaussian())

ols_tiles <- fixest::feglm(avg_daily_tiles ~ d_log_fwc_100k | as.factor(userid) +
                             fw_median_hh_inc_grp_X_month + fw_pop_density_grp_X_month + fw_frac_urban_grp_X_month +
                             zcta_interaction_group_X_month,
                           data=dat, fixef.rm = "both", family = gaussian())

poisson_tiles <- fixest::feglm(avg_daily_tiles ~ d_log_fwc_100k | as.factor(userid) +
                                 fw_median_hh_inc_grp_X_month + fw_pop_density_grp_X_month + fw_frac_urban_grp_X_month +
                                 zcta_interaction_group_X_month,
                               data=dat, fixef.rm = "both", family = "poisson")

# Make a list of the models
models <- list(ols_sanity, ols_tiles_diff_no_userid, ols_tiles_diff, ols_tiles, poisson_tiles)

# Add the y-means
model_y_mean <- function(model){ mean(model$residuals + model$fitted.values) }
y_means <- sapply(models, model_y_mean)

out <- etable(models,
            cluster = ~home_zcta,
            headers = c("OLS Sanity",
                        "OLS Tiles Diff No Userid", "OLS Tiles Diff",
                        "OLS Tiles", "Poisson Tiles"),
            extralines=list("_Y-Mean"=y_means))

names(out)[1] <- "_"

out <- out %>% as_tibble()

write_csv(out, "output/tables/fixest_robustness_regs.csv")
